<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">


<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
    
    <title>Raster Layers &mdash; Python GDAL/OGR Cookbook 1.0 documentation</title>
    
    <link rel="stylesheet" href="_static/classic.css" type="text/css" />
    <link rel="stylesheet" href="_static/pygments.css" type="text/css" />
    
    <script type="text/javascript">
      var DOCUMENTATION_OPTIONS = {
        URL_ROOT:    './',
        VERSION:     '1.0',
        COLLAPSE_INDEX: false,
        FILE_SUFFIX: '.html',
        HAS_SOURCE:  true
      };
    </script>
    <script type="text/javascript" src="_static/jquery.js"></script>
    <script type="text/javascript" src="_static/underscore.js"></script>
    <script type="text/javascript" src="_static/doctools.js"></script>
    <link rel="top" title="Python GDAL/OGR Cookbook 1.0 documentation" href="index.html" />
    <link rel="next" title="Projection" href="projection.html" />
    <link rel="prev" title="Vector Layers" href="vector_layers.html" /> 
  </head>
  <body role="document">
    <div class="related" role="navigation" aria-label="related navigation">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="genindex.html" title="General Index"
             accesskey="I">index</a></li>
        <li class="right" >
          <a href="projection.html" title="Projection"
             accesskey="N">next</a> |</li>
        <li class="right" >
          <a href="vector_layers.html" title="Vector Layers"
             accesskey="P">previous</a> |</li>
        <li class="nav-item nav-item-0"><a href="index.html">Python GDAL/OGR Cookbook 1.0 documentation</a> &raquo;</li> 
      </ul>
    </div>  

    <div class="document">
      <div class="documentwrapper">
        <div class="bodywrapper">
          <div class="body" role="main">
            
  <div class="section" id="raster-layers">
<h1>Raster Layers<a class="headerlink" href="#raster-layers" title="Permalink to this headline">¶</a></h1>
<div class="section" id="close-a-raster-dataset">
<h2>Close a raster dataset<a class="headerlink" href="#close-a-raster-dataset" title="Permalink to this headline">¶</a></h2>
<p>This recipe shows how to close a raster dataset.
It is useful in the middle of a script, to recover the resources held by accessing the dataset, remove file locks, etc. It is not necessary at the end of the script, as the Python garbage collector will do the same thing automatically when the script exits.
Further discussion of this topic can be found in this <a class="reference external" href="http://gis.stackexchange.com/questions/80366/why-close-a-dataset-in-gdal-python">GIS StackExchange question</a>.</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">import</span> <span class="nn">gdal</span>

<span class="c"># open dataset</span>
<span class="n">ds</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span><span class="s">&#39;test.tif&#39;</span><span class="p">)</span>

<span class="c"># close dataset</span>
<span class="n">ds</span> <span class="o">=</span> <span class="bp">None</span>
</pre></div>
</div>
</div>
<div class="section" id="get-raster-metadata">
<h2>Get Raster Metadata<a class="headerlink" href="#get-raster-metadata" title="Permalink to this headline">¶</a></h2>
<p>Get raster metadata for quick-and-dirty resolution checks</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">from</span> <span class="nn">osgeo</span> <span class="kn">import</span> <span class="n">gdal</span>
<span class="n">gtif</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span> <span class="s">&quot;INPUT.tif&quot;</span> <span class="p">)</span>
<span class="k">print</span> <span class="n">gtif</span><span class="o">.</span><span class="n">GetMetadata</span><span class="p">()</span>
</pre></div>
</div>
</div>
<div class="section" id="get-raster-band">
<h2>Get Raster Band<a class="headerlink" href="#get-raster-band" title="Permalink to this headline">¶</a></h2>
<p>Get a raster band. Notice how we are handling runtime errors this function might throw.</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">from</span> <span class="nn">osgeo</span> <span class="kn">import</span> <span class="n">gdal</span>
<span class="kn">import</span> <span class="nn">sys</span>
<span class="c"># this allows GDAL to throw Python Exceptions</span>
<span class="n">gdal</span><span class="o">.</span><span class="n">UseExceptions</span><span class="p">()</span>

<span class="k">try</span><span class="p">:</span>
    <span class="n">src_ds</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span> <span class="s">&quot;INPUT.tif&quot;</span> <span class="p">)</span>
<span class="k">except</span> <span class="ne">RuntimeError</span><span class="p">,</span> <span class="n">e</span><span class="p">:</span>
    <span class="k">print</span> <span class="s">&#39;Unable to open INPUT.tif&#39;</span>
    <span class="k">print</span> <span class="n">e</span>
    <span class="n">sys</span><span class="o">.</span><span class="n">exit</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>

<span class="k">try</span><span class="p">:</span>
    <span class="n">srcband</span> <span class="o">=</span> <span class="n">src_ds</span><span class="o">.</span><span class="n">GetRasterBand</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>
<span class="k">except</span> <span class="ne">RuntimeError</span><span class="p">,</span> <span class="n">e</span><span class="p">:</span>
    <span class="c"># for example, try GetRasterBand(10)</span>
    <span class="k">print</span> <span class="s">&#39;Band ( </span><span class="si">%i</span><span class="s"> ) not found&#39;</span> <span class="o">%</span> <span class="n">band_num</span>
    <span class="k">print</span> <span class="n">e</span>
    <span class="n">sys</span><span class="o">.</span><span class="n">exit</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="section" id="loop-through-all-raster-bands">
<h2>Loop Through All Raster Bands<a class="headerlink" href="#loop-through-all-raster-bands" title="Permalink to this headline">¶</a></h2>
<p>Loop through all raster bands and do something useful like listing band statistics.</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">from</span> <span class="nn">osgeo</span> <span class="kn">import</span> <span class="n">gdal</span>
<span class="kn">import</span> <span class="nn">sys</span>

<span class="n">src_ds</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span> <span class="s">&quot;INPUT.tif&quot;</span> <span class="p">)</span>
<span class="k">if</span> <span class="n">src_ds</span> <span class="ow">is</span> <span class="bp">None</span><span class="p">:</span>
    <span class="k">print</span> <span class="s">&#39;Unable to open INPUT.tif&#39;</span>
    <span class="n">sys</span><span class="o">.</span><span class="n">exit</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>

<span class="k">print</span> <span class="s">&quot;[ RASTER BAND COUNT ]: &quot;</span><span class="p">,</span> <span class="n">src_ds</span><span class="o">.</span><span class="n">RasterCount</span>
<span class="k">for</span> <span class="n">band</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span> <span class="n">src_ds</span><span class="o">.</span><span class="n">RasterCount</span> <span class="p">):</span>
    <span class="n">band</span> <span class="o">+=</span> <span class="mi">1</span>
    <span class="k">print</span> <span class="s">&quot;[ GETTING BAND ]: &quot;</span><span class="p">,</span> <span class="n">band</span>
    <span class="n">srcband</span> <span class="o">=</span> <span class="n">src_ds</span><span class="o">.</span><span class="n">GetRasterBand</span><span class="p">(</span><span class="n">band</span><span class="p">)</span>
    <span class="k">if</span> <span class="n">srcband</span> <span class="ow">is</span> <span class="bp">None</span><span class="p">:</span>
        <span class="k">continue</span>

    <span class="n">stats</span> <span class="o">=</span> <span class="n">srcband</span><span class="o">.</span><span class="n">GetStatistics</span><span class="p">(</span> <span class="bp">True</span><span class="p">,</span> <span class="bp">True</span> <span class="p">)</span>
    <span class="k">if</span> <span class="n">stats</span> <span class="ow">is</span> <span class="bp">None</span><span class="p">:</span>
        <span class="k">continue</span>

    <span class="k">print</span> <span class="s">&quot;[ STATS ] =  Minimum=</span><span class="si">%.3f</span><span class="s">, Maximum=</span><span class="si">%.3f</span><span class="s">, Mean=</span><span class="si">%.3f</span><span class="s">, StdDev=</span><span class="si">%.3f</span><span class="s">&quot;</span> <span class="o">%</span> <span class="p">(</span> \
                <span class="n">stats</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span> <span class="n">stats</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span> <span class="n">stats</span><span class="p">[</span><span class="mi">2</span><span class="p">],</span> <span class="n">stats</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> <span class="p">)</span>
</pre></div>
</div>
</div>
<div class="section" id="get-raster-band-information">
<h2>Get Raster Band Information<a class="headerlink" href="#get-raster-band-information" title="Permalink to this headline">¶</a></h2>
<p>Using the documentation on the <a class="reference external" href="http://gdal.org/python/osgeo.gdal.Band-class.html">Band API</a> we can
write a script that dumps out single band information</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">from</span> <span class="nn">osgeo</span> <span class="kn">import</span> <span class="n">gdal</span>
<span class="kn">import</span> <span class="nn">sys</span>
<span class="n">gdal</span><span class="o">.</span><span class="n">UseExceptions</span><span class="p">()</span>

<span class="k">def</span> <span class="nf">Usage</span><span class="p">():</span>
    <span class="k">print</span><span class="p">(</span><span class="s">&quot;&quot;&quot;</span>
<span class="s">    $ getrasterband.py [ band number ] input-raster</span>
<span class="s">    &quot;&quot;&quot;</span><span class="p">)</span>
    <span class="n">sys</span><span class="o">.</span><span class="n">exit</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>

<span class="k">def</span> <span class="nf">main</span><span class="p">(</span> <span class="n">band_num</span><span class="p">,</span> <span class="n">input_file</span> <span class="p">):</span>
    <span class="n">src_ds</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span> <span class="n">input_file</span> <span class="p">)</span>
    <span class="k">if</span> <span class="n">src_ds</span> <span class="ow">is</span> <span class="bp">None</span><span class="p">:</span>
        <span class="k">print</span> <span class="s">&#39;Unable to open </span><span class="si">%s</span><span class="s">&#39;</span> <span class="o">%</span> <span class="n">input_file</span>
        <span class="n">sys</span><span class="o">.</span><span class="n">exit</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>

    <span class="k">try</span><span class="p">:</span>
        <span class="n">srcband</span> <span class="o">=</span> <span class="n">src_ds</span><span class="o">.</span><span class="n">GetRasterBand</span><span class="p">(</span><span class="n">band_num</span><span class="p">)</span>
    <span class="k">except</span> <span class="ne">RuntimeError</span><span class="p">,</span> <span class="n">e</span><span class="p">:</span>
        <span class="k">print</span> <span class="s">&#39;No band </span><span class="si">%i</span><span class="s"> found&#39;</span> <span class="o">%</span> <span class="n">band_num</span>
        <span class="k">print</span> <span class="n">e</span>
        <span class="n">sys</span><span class="o">.</span><span class="n">exit</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>


    <span class="k">print</span> <span class="s">&quot;[ NO DATA VALUE ] = &quot;</span><span class="p">,</span> <span class="n">srcband</span><span class="o">.</span><span class="n">GetNoDataValue</span><span class="p">()</span>
    <span class="k">print</span> <span class="s">&quot;[ MIN ] = &quot;</span><span class="p">,</span> <span class="n">srcband</span><span class="o">.</span><span class="n">GetMinimum</span><span class="p">()</span>
    <span class="k">print</span> <span class="s">&quot;[ MAX ] = &quot;</span><span class="p">,</span> <span class="n">srcband</span><span class="o">.</span><span class="n">GetMaximum</span><span class="p">()</span>
    <span class="k">print</span> <span class="s">&quot;[ SCALE ] = &quot;</span><span class="p">,</span> <span class="n">srcband</span><span class="o">.</span><span class="n">GetScale</span><span class="p">()</span>
    <span class="k">print</span> <span class="s">&quot;[ UNIT TYPE ] = &quot;</span><span class="p">,</span> <span class="n">srcband</span><span class="o">.</span><span class="n">GetUnitType</span><span class="p">()</span>
    <span class="n">ctable</span> <span class="o">=</span> <span class="n">srcband</span><span class="o">.</span><span class="n">GetColorTable</span><span class="p">()</span>

    <span class="k">if</span> <span class="n">ctable</span> <span class="ow">is</span> <span class="bp">None</span><span class="p">:</span>
        <span class="k">print</span> <span class="s">&#39;No ColorTable found&#39;</span>
        <span class="n">sys</span><span class="o">.</span><span class="n">exit</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>

    <span class="k">print</span> <span class="s">&quot;[ COLOR TABLE COUNT ] = &quot;</span><span class="p">,</span> <span class="n">ctable</span><span class="o">.</span><span class="n">GetCount</span><span class="p">()</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span> <span class="mi">0</span><span class="p">,</span> <span class="n">ctable</span><span class="o">.</span><span class="n">GetCount</span><span class="p">()</span> <span class="p">):</span>
        <span class="n">entry</span> <span class="o">=</span> <span class="n">ctable</span><span class="o">.</span><span class="n">GetColorEntry</span><span class="p">(</span> <span class="n">i</span> <span class="p">)</span>
        <span class="k">if</span> <span class="ow">not</span> <span class="n">entry</span><span class="p">:</span>
            <span class="k">continue</span>
        <span class="k">print</span> <span class="s">&quot;[ COLOR ENTRY RGB ] = &quot;</span><span class="p">,</span> <span class="n">ctable</span><span class="o">.</span><span class="n">GetColorEntryAsRGB</span><span class="p">(</span> <span class="n">i</span><span class="p">,</span> <span class="n">entry</span> <span class="p">)</span>

<span class="k">if</span> <span class="n">__name__</span> <span class="o">==</span> <span class="s">&#39;__main__&#39;</span><span class="p">:</span>

    <span class="k">if</span> <span class="nb">len</span><span class="p">(</span> <span class="n">sys</span><span class="o">.</span><span class="n">argv</span> <span class="p">)</span> <span class="o">&lt;</span> <span class="mi">3</span><span class="p">:</span>
        <span class="k">print</span> <span class="s">&quot;&quot;&quot;</span>
<span class="s">        [ ERROR ] you must supply at least two arguments:</span>
<span class="s">        1) the band number to retrieve and 2) input raster</span>
<span class="s">        &quot;&quot;&quot;</span>
        <span class="n">Usage</span><span class="p">()</span>

    <span class="n">main</span><span class="p">(</span> <span class="nb">int</span><span class="p">(</span><span class="n">sys</span><span class="o">.</span><span class="n">argv</span><span class="p">[</span><span class="mi">1</span><span class="p">]),</span> <span class="n">sys</span><span class="o">.</span><span class="n">argv</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span> <span class="p">)</span>
</pre></div>
</div>
</div>
<div class="section" id="polygonize-a-raster-band">
<h2>Polygonize a Raster Band<a class="headerlink" href="#polygonize-a-raster-band" title="Permalink to this headline">¶</a></h2>
<p>Turn a single raster band into a vector polygon!</p>
<p>If you haven&#8217;t before, notice that some of the <a class="reference external" href="http://www.gdal.org/gdal_utilities.html">gdal utilties</a>
are actually Python scripts. Go find them on your computer, read the source code and mine them for API tricks.
It turns out the <a class="reference external" href="http://www.gdal.org/gdal_polygonize.html">gdal_polygonize utility</a>
just wraps a call to <a class="reference external" href="http://www.gdal.org/gdal__alg_8h.html#a3f522a9035d3512b5d414fb4752671b1">GDALFPolygonize</a>
so writing your own hacky polygonize Python script is pretty easy.</p>
<p>The raster we are going to polygonize:</p>
<img alt="_images/input_tif.jpg" src="_images/input_tif.jpg" />
<div class="highlight-python"><div class="highlight"><pre><span class="kn">from</span> <span class="nn">osgeo</span> <span class="kn">import</span> <span class="n">gdal</span><span class="p">,</span> <span class="n">ogr</span>
<span class="kn">import</span> <span class="nn">sys</span>
<span class="c"># this allows GDAL to throw Python Exceptions</span>
<span class="n">gdal</span><span class="o">.</span><span class="n">UseExceptions</span><span class="p">()</span>

<span class="c">#</span>
<span class="c">#  get raster datasource</span>
<span class="c">#</span>
<span class="n">src_ds</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span> <span class="s">&quot;INPUT.tif&quot;</span> <span class="p">)</span>
<span class="k">if</span> <span class="n">src_ds</span> <span class="ow">is</span> <span class="bp">None</span><span class="p">:</span>
    <span class="k">print</span> <span class="s">&#39;Unable to open </span><span class="si">%s</span><span class="s">&#39;</span> <span class="o">%</span> <span class="n">src_filename</span>
    <span class="n">sys</span><span class="o">.</span><span class="n">exit</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>

<span class="k">try</span><span class="p">:</span>
    <span class="n">srcband</span> <span class="o">=</span> <span class="n">src_ds</span><span class="o">.</span><span class="n">GetRasterBand</span><span class="p">(</span><span class="mi">3</span><span class="p">)</span>
<span class="k">except</span> <span class="ne">RuntimeError</span><span class="p">,</span> <span class="n">e</span><span class="p">:</span>
    <span class="c"># for example, try GetRasterBand(10)</span>
    <span class="k">print</span> <span class="s">&#39;Band ( </span><span class="si">%i</span><span class="s"> ) not found&#39;</span> <span class="o">%</span> <span class="n">band_num</span>
    <span class="k">print</span> <span class="n">e</span>
    <span class="n">sys</span><span class="o">.</span><span class="n">exit</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>

<span class="c">#</span>
<span class="c">#  create output datasource</span>
<span class="c">#</span>
<span class="n">dst_layername</span> <span class="o">=</span> <span class="s">&quot;POLYGONIZED_STUFF&quot;</span>
<span class="n">drv</span> <span class="o">=</span> <span class="n">ogr</span><span class="o">.</span><span class="n">GetDriverByName</span><span class="p">(</span><span class="s">&quot;ESRI Shapefile&quot;</span><span class="p">)</span>
<span class="n">dst_ds</span> <span class="o">=</span> <span class="n">drv</span><span class="o">.</span><span class="n">CreateDataSource</span><span class="p">(</span> <span class="n">dst_layername</span> <span class="o">+</span> <span class="s">&quot;.shp&quot;</span> <span class="p">)</span>
<span class="n">dst_layer</span> <span class="o">=</span> <span class="n">dst_ds</span><span class="o">.</span><span class="n">CreateLayer</span><span class="p">(</span><span class="n">dst_layername</span><span class="p">,</span> <span class="n">srs</span> <span class="o">=</span> <span class="bp">None</span> <span class="p">)</span>

<span class="n">gdal</span><span class="o">.</span><span class="n">Polygonize</span><span class="p">(</span> <span class="n">srcband</span><span class="p">,</span> <span class="bp">None</span><span class="p">,</span> <span class="n">dst_layer</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="p">[],</span> <span class="n">callback</span><span class="o">=</span><span class="bp">None</span> <span class="p">)</span>
</pre></div>
</div>
<p>A closeup of the vector output &#8211; this is northern Italy:</p>
<img alt="_images/polygonize_band3.png" src="_images/polygonize_band3.png" />
</div>
<div class="section" id="convert-an-ogr-file-to-a-raster">
<h2>Convert an OGR File to a Raster<a class="headerlink" href="#convert-an-ogr-file-to-a-raster" title="Permalink to this headline">¶</a></h2>
<p>This recipe takes in a OGR file (e.g. shapefile) and creates a new raster Tiff file based on the shapefile.</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">from</span> <span class="nn">osgeo</span> <span class="kn">import</span> <span class="n">gdal</span><span class="p">,</span> <span class="n">ogr</span>

<span class="c"># Define pixel_size and NoData value of new raster</span>
<span class="n">pixel_size</span> <span class="o">=</span> <span class="mi">25</span>
<span class="n">NoData_value</span> <span class="o">=</span> <span class="o">-</span><span class="mi">9999</span>

<span class="c"># Filename of input OGR file</span>
<span class="n">vector_fn</span> <span class="o">=</span> <span class="s">&#39;test.shp&#39;</span>

<span class="c"># Filename of the raster Tiff that will be created</span>
<span class="n">raster_fn</span> <span class="o">=</span> <span class="s">&#39;test.tif&#39;</span>

<span class="c"># Open the data source and read in the extent</span>
<span class="n">source_ds</span> <span class="o">=</span> <span class="n">ogr</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span><span class="n">vector_fn</span><span class="p">)</span>
<span class="n">source_layer</span> <span class="o">=</span> <span class="n">source_ds</span><span class="o">.</span><span class="n">GetLayer</span><span class="p">()</span>
<span class="n">x_min</span><span class="p">,</span> <span class="n">x_max</span><span class="p">,</span> <span class="n">y_min</span><span class="p">,</span> <span class="n">y_max</span> <span class="o">=</span> <span class="n">source_layer</span><span class="o">.</span><span class="n">GetExtent</span><span class="p">()</span>

<span class="c"># Create the destination data source</span>
<span class="n">x_res</span> <span class="o">=</span> <span class="nb">int</span><span class="p">((</span><span class="n">x_max</span> <span class="o">-</span> <span class="n">x_min</span><span class="p">)</span> <span class="o">/</span> <span class="n">pixel_size</span><span class="p">)</span>
<span class="n">y_res</span> <span class="o">=</span> <span class="nb">int</span><span class="p">((</span><span class="n">y_max</span> <span class="o">-</span> <span class="n">y_min</span><span class="p">)</span> <span class="o">/</span> <span class="n">pixel_size</span><span class="p">)</span>
<span class="n">target_ds</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">GetDriverByName</span><span class="p">(</span><span class="s">&#39;GTiff&#39;</span><span class="p">)</span><span class="o">.</span><span class="n">Create</span><span class="p">(</span><span class="n">raster_fn</span><span class="p">,</span> <span class="n">x_res</span><span class="p">,</span> <span class="n">y_res</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="n">gdal</span><span class="o">.</span><span class="n">GDT_Byte</span><span class="p">)</span>
<span class="n">target_ds</span><span class="o">.</span><span class="n">SetGeoTransform</span><span class="p">((</span><span class="n">x_min</span><span class="p">,</span> <span class="n">pixel_size</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="n">y_max</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="o">-</span><span class="n">pixel_size</span><span class="p">))</span>
<span class="n">band</span> <span class="o">=</span> <span class="n">target_ds</span><span class="o">.</span><span class="n">GetRasterBand</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>
<span class="n">band</span><span class="o">.</span><span class="n">SetNoDataValue</span><span class="p">(</span><span class="n">NoData_value</span><span class="p">)</span>

<span class="c"># Rasterize</span>
<span class="n">gdal</span><span class="o">.</span><span class="n">RasterizeLayer</span><span class="p">(</span><span class="n">target_ds</span><span class="p">,</span> <span class="p">[</span><span class="mi">1</span><span class="p">],</span> <span class="n">source_layer</span><span class="p">,</span> <span class="n">burn_values</span><span class="o">=</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span>
</pre></div>
</div>
</div>
<div class="section" id="clip-a-geotiff-with-shapefile">
<h2>Clip a GeoTiff with Shapefile<a class="headerlink" href="#clip-a-geotiff-with-shapefile" title="Permalink to this headline">¶</a></h2>
<p>Let&#8217;s use some <a class="reference external" href="http://www.naturalearthdata.com/downloads">Natural Earth data</a> and clip a <a class="reference external" href="http://www.naturalearthdata.com/downloads/10m-cross-blend-hypso/cross-blended-hypso-with-relief-water-drains-and-ocean-bottom/">10m relief geotiff</a> with the <a class="reference external" href="http://www.naturalearthdata.com/downloads/10m-cultural-vectors/timezones/">Europe/Paris timezone polygon</a>. Most of the following workflow came from this <a class="reference external" href="http://geospatialpython.com/2011/02/clip-raster-using-shapefile.html">geospatialpython post</a> . However, the source code on that site assumes your clipping polygon <strong>is</strong> the same extent as the input geotiff. If it is not, then your clipped geotiff will take the input geotiff&#8217;s extent, which will be incorrect. The modified script below takes this into account and sets the correct x,y offsets for the clipped geotiff. Note, in the following example we are assuming you have the <a class="reference external" href="http://www.pythonware.com/products/pil/">Python Imaging Library</a> installed.</p>
<p>Before Image: the input Natural Earth 10m geotiff with the timezone overlay we want to clip out:</p>
<img alt="_images/clip_raster_before.png" src="_images/clip_raster_before.png" />
<div class="highlight-bash"><div class="highlight"><pre>from osgeo import gdal, gdalnumeric, ogr, osr
import Image, ImageDraw
import os, sys
gdal.UseExceptions<span class="o">()</span>


<span class="c"># This function will convert the rasterized clipper shapefile</span>
<span class="c"># to a mask for use within GDAL.</span>
def imageToArray<span class="o">(</span>i<span class="o">)</span>:
    <span class="s2">&quot;&quot;&quot;</span>
<span class="s2">    Converts a Python Imaging Library array to a</span>
<span class="s2">    gdalnumeric image.</span>
<span class="s2">    &quot;&quot;&quot;</span>
    <span class="nv">a</span><span class="o">=</span>gdalnumeric.fromstring<span class="o">(</span>i.tostring<span class="o">()</span>,<span class="s1">&#39;b&#39;</span><span class="o">)</span>
    a.shape<span class="o">=</span>i.im.size<span class="o">[</span>1<span class="o">]</span>, i.im.size<span class="o">[</span>0<span class="o">]</span>
    <span class="k">return</span> a

def arrayToImage<span class="o">(</span>a<span class="o">)</span>:
    <span class="s2">&quot;&quot;&quot;</span>
<span class="s2">    Converts a gdalnumeric array to a</span>
<span class="s2">    Python Imaging Library Image.</span>
<span class="s2">    &quot;&quot;&quot;</span>
    <span class="nv">i</span><span class="o">=</span>Image.fromstring<span class="o">(</span><span class="s1">&#39;L&#39;</span>,<span class="o">(</span>a.shape<span class="o">[</span>1<span class="o">]</span>,a.shape<span class="o">[</span>0<span class="o">])</span>,
            <span class="o">(</span>a.astype<span class="o">(</span><span class="s1">&#39;b&#39;</span><span class="o">))</span>.tostring<span class="o">())</span>
    <span class="k">return</span> i

def world2Pixel<span class="o">(</span>geoMatrix, x, y<span class="o">)</span>:
  <span class="s2">&quot;&quot;&quot;</span>
<span class="s2">  Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate</span>
<span class="s2">  the pixel location of a geospatial coordinate</span>
<span class="s2">  &quot;&quot;&quot;</span>
  <span class="nv">ulX</span> <span class="o">=</span> geoMatrix<span class="o">[</span>0<span class="o">]</span>
  <span class="nv">ulY</span> <span class="o">=</span> geoMatrix<span class="o">[</span>3<span class="o">]</span>
  <span class="nv">xDist</span> <span class="o">=</span> geoMatrix<span class="o">[</span>1<span class="o">]</span>
  <span class="nv">yDist</span> <span class="o">=</span> geoMatrix<span class="o">[</span>5<span class="o">]</span>
  <span class="nv">rtnX</span> <span class="o">=</span> geoMatrix<span class="o">[</span>2<span class="o">]</span>
  <span class="nv">rtnY</span> <span class="o">=</span> geoMatrix<span class="o">[</span>4<span class="o">]</span>
  <span class="nv">pixel</span> <span class="o">=</span> int<span class="o">((</span>x - ulX<span class="o">)</span> / xDist<span class="o">)</span>
  <span class="nv">line</span> <span class="o">=</span> int<span class="o">((</span>ulY - y<span class="o">)</span> / xDist<span class="o">)</span>
  <span class="k">return</span> <span class="o">(</span>pixel, line<span class="o">)</span>

<span class="c">#</span>
<span class="c">#  EDIT: this is basically an overloaded</span>
<span class="c">#  version of the gdal_array.OpenArray passing in xoff, yoff explicitly</span>
<span class="c">#  so we can pass these params off to CopyDatasetInfo</span>
<span class="c">#</span>
def OpenArray<span class="o">(</span> array, <span class="nv">prototype_ds</span> <span class="o">=</span> None, <span class="nv">xoff</span><span class="o">=</span>0, <span class="nv">yoff</span><span class="o">=</span><span class="m">0</span> <span class="o">)</span>:
    <span class="nv">ds</span> <span class="o">=</span> gdal.Open<span class="o">(</span> gdalnumeric.GetArrayFilename<span class="o">(</span>array<span class="o">)</span> <span class="o">)</span>

    <span class="k">if</span> ds is not None and prototype_ds is not None:
        <span class="k">if</span> <span class="nb">type</span><span class="o">(</span>prototype_ds<span class="o">)</span>.__name__ <span class="o">==</span> <span class="s1">&#39;str&#39;</span>:
            <span class="nv">prototype_ds</span> <span class="o">=</span> gdal.Open<span class="o">(</span> prototype_ds <span class="o">)</span>
        <span class="k">if</span> prototype_ds is not None:
            gdalnumeric.CopyDatasetInfo<span class="o">(</span> prototype_ds, ds, <span class="nv">xoff</span><span class="o">=</span>xoff, <span class="nv">yoff</span><span class="o">=</span>yoff <span class="o">)</span>
    <span class="k">return</span> ds

def histogram<span class="o">(</span>a, <span class="nv">bins</span><span class="o">=</span>range<span class="o">(</span>0,256<span class="o">))</span>:
  <span class="s2">&quot;&quot;&quot;</span>
<span class="s2">  Histogram function for multi-dimensional array.</span>
<span class="s2">  a = array</span>
<span class="s2">  bins = range of numbers to match</span>
<span class="s2">  &quot;&quot;&quot;</span>
  <span class="nv">fa</span> <span class="o">=</span> a.flat
  <span class="nv">n</span> <span class="o">=</span> gdalnumeric.searchsorted<span class="o">(</span>gdalnumeric.sort<span class="o">(</span>fa<span class="o">)</span>, bins<span class="o">)</span>
  <span class="nv">n</span> <span class="o">=</span> gdalnumeric.concatenate<span class="o">([</span>n, <span class="o">[</span>len<span class="o">(</span>fa<span class="o">)]])</span>
  <span class="nv">hist</span> <span class="o">=</span> n<span class="o">[</span>1:<span class="o">]</span>-n<span class="o">[</span>:-1<span class="o">]</span>
  <span class="k">return</span> hist

def stretch<span class="o">(</span>a<span class="o">)</span>:
  <span class="s2">&quot;&quot;&quot;</span>
<span class="s2">  Performs a histogram stretch on a gdalnumeric array image.</span>
<span class="s2">  &quot;&quot;&quot;</span>
  <span class="nv">hist</span> <span class="o">=</span> histogram<span class="o">(</span>a<span class="o">)</span>
  <span class="nv">im</span> <span class="o">=</span> arrayToImage<span class="o">(</span>a<span class="o">)</span>
  <span class="nv">lut</span> <span class="o">=</span> <span class="o">[]</span>
  <span class="k">for</span> b in range<span class="o">(</span>0, len<span class="o">(</span>hist<span class="o">)</span>, 256<span class="o">)</span>:
    <span class="c"># step size</span>
    <span class="nv">step</span> <span class="o">=</span> reduce<span class="o">(</span>operator.add, hist<span class="o">[</span>b:b+256<span class="o">])</span> / 255
    <span class="c"># create equalization lookup table</span>
    <span class="nv">n</span> <span class="o">=</span> 0
    <span class="k">for</span> i in range<span class="o">(</span>256<span class="o">)</span>:
      lut.append<span class="o">(</span>n / step<span class="o">)</span>
      <span class="nv">n</span> <span class="o">=</span> n + hist<span class="o">[</span>i+b<span class="o">]</span>
  <span class="nv">im</span> <span class="o">=</span> im.point<span class="o">(</span>lut<span class="o">)</span>
  <span class="k">return</span> imageToArray<span class="o">(</span>im<span class="o">)</span>

def main<span class="o">(</span> shapefile_path, raster_path <span class="o">)</span>:
    <span class="c"># Load the source data as a gdalnumeric array</span>
    <span class="nv">srcArray</span> <span class="o">=</span> gdalnumeric.LoadFile<span class="o">(</span>raster_path<span class="o">)</span>

    <span class="c"># Also load as a gdal image to get geotransform</span>
    <span class="c"># (world file) info</span>
    <span class="nv">srcImage</span> <span class="o">=</span> gdal.Open<span class="o">(</span>raster_path<span class="o">)</span>
    <span class="nv">geoTrans</span> <span class="o">=</span> srcImage.GetGeoTransform<span class="o">()</span>

    <span class="c"># Create an OGR layer from a boundary shapefile</span>
    <span class="nv">shapef</span> <span class="o">=</span> ogr.Open<span class="o">(</span>shapefile_path<span class="o">)</span>
    <span class="nv">lyr</span> <span class="o">=</span> shapef.GetLayer<span class="o">(</span> os.path.split<span class="o">(</span> os.path.splitext<span class="o">(</span> shapefile_path <span class="o">)[</span>0<span class="o">]</span> <span class="o">)[</span>1<span class="o">]</span> <span class="o">)</span>
    <span class="nv">poly</span> <span class="o">=</span> lyr.GetNextFeature<span class="o">()</span>

    <span class="c"># Convert the layer extent to image pixel coordinates</span>
    minX, maxX, minY, <span class="nv">maxY</span> <span class="o">=</span> lyr.GetExtent<span class="o">()</span>
    ulX, <span class="nv">ulY</span> <span class="o">=</span> world2Pixel<span class="o">(</span>geoTrans, minX, maxY<span class="o">)</span>
    lrX, <span class="nv">lrY</span> <span class="o">=</span> world2Pixel<span class="o">(</span>geoTrans, maxX, minY<span class="o">)</span>

    <span class="c"># Calculate the pixel size of the new image</span>
    <span class="nv">pxWidth</span> <span class="o">=</span> int<span class="o">(</span>lrX - ulX<span class="o">)</span>
    <span class="nv">pxHeight</span> <span class="o">=</span> int<span class="o">(</span>lrY - ulY<span class="o">)</span>

    <span class="nv">clip</span> <span class="o">=</span> srcArray<span class="o">[</span>:, ulY:lrY, ulX:lrX<span class="o">]</span>

    <span class="c">#</span>
    <span class="c"># EDIT: create pixel offset to pass to new image Projection info</span>
    <span class="c">#</span>
    <span class="nv">xoffset</span> <span class="o">=</span>  ulX
    <span class="nv">yoffset</span> <span class="o">=</span>  ulY
    print <span class="s2">&quot;Xoffset, Yoffset = ( %f, %f )&quot;</span> % <span class="o">(</span> xoffset, yoffset <span class="o">)</span>

    <span class="c"># Create a new geomatrix for the image</span>
    <span class="nv">geoTrans</span> <span class="o">=</span> list<span class="o">(</span>geoTrans<span class="o">)</span>
    geoTrans<span class="o">[</span>0<span class="o">]</span> <span class="o">=</span> minX
    geoTrans<span class="o">[</span>3<span class="o">]</span> <span class="o">=</span> maxY

    <span class="c"># Map points to pixels for drawing the</span>
    <span class="c"># boundary on a blank 8-bit,</span>
    <span class="c"># black and white, mask image.</span>
    <span class="nv">points</span> <span class="o">=</span> <span class="o">[]</span>
    <span class="nv">pixels</span> <span class="o">=</span> <span class="o">[]</span>
    <span class="nv">geom</span> <span class="o">=</span> poly.GetGeometryRef<span class="o">()</span>
    <span class="nv">pts</span> <span class="o">=</span> geom.GetGeometryRef<span class="o">(</span>0<span class="o">)</span>
    <span class="k">for</span> p in range<span class="o">(</span>pts.GetPointCount<span class="o">())</span>:
      points.append<span class="o">((</span>pts.GetX<span class="o">(</span>p<span class="o">)</span>, pts.GetY<span class="o">(</span>p<span class="o">)))</span>
    <span class="k">for</span> p in points:
      pixels.append<span class="o">(</span>world2Pixel<span class="o">(</span>geoTrans, p<span class="o">[</span>0<span class="o">]</span>, p<span class="o">[</span>1<span class="o">]))</span>
    <span class="nv">rasterPoly</span> <span class="o">=</span> Image.new<span class="o">(</span><span class="s2">&quot;L&quot;</span>, <span class="o">(</span>pxWidth, pxHeight<span class="o">)</span>, 1<span class="o">)</span>
    <span class="nv">rasterize</span> <span class="o">=</span> ImageDraw.Draw<span class="o">(</span>rasterPoly<span class="o">)</span>
    rasterize.polygon<span class="o">(</span>pixels, 0<span class="o">)</span>
    <span class="nv">mask</span> <span class="o">=</span> imageToArray<span class="o">(</span>rasterPoly<span class="o">)</span>

    <span class="c"># Clip the image using the mask</span>
    <span class="nv">clip</span> <span class="o">=</span> gdalnumeric.choose<span class="o">(</span>mask, <span class="se">\</span>
        <span class="o">(</span>clip, 0<span class="o">))</span>.astype<span class="o">(</span>gdalnumeric.uint8<span class="o">)</span>

    <span class="c"># This image has 3 bands so we stretch each one to make them</span>
    <span class="c"># visually brighter</span>
    <span class="k">for</span> i in range<span class="o">(</span>3<span class="o">)</span>:
      clip<span class="o">[</span>i,:,:<span class="o">]</span> <span class="o">=</span> stretch<span class="o">(</span>clip<span class="o">[</span>i,:,:<span class="o">])</span>

    <span class="c"># Save new tiff</span>
    <span class="c">#</span>
    <span class="c">#  EDIT: instead of SaveArray, let&#39;s break all the</span>
    <span class="c">#  SaveArray steps out more explicity so</span>
    <span class="c">#  we can overwrite the offset of the destination</span>
    <span class="c">#  raster</span>
    <span class="c">#</span>
    <span class="c">### the old way using SaveArray</span>
    <span class="c">#</span>
    <span class="c"># gdalnumeric.SaveArray(clip, &quot;OUTPUT.tif&quot;, format=&quot;GTiff&quot;, prototype=raster_path)</span>
    <span class="c">#</span>
    <span class="c">###</span>
    <span class="c">#</span>
    <span class="nv">gtiffDriver</span> <span class="o">=</span> gdal.GetDriverByName<span class="o">(</span> <span class="s1">&#39;GTiff&#39;</span> <span class="o">)</span>
    <span class="k">if</span> gtiffDriver is None:
        raise ValueError<span class="o">(</span><span class="s2">&quot;Can&#39;t find GeoTiff Driver&quot;</span><span class="o">)</span>
    gtiffDriver.CreateCopy<span class="o">(</span> <span class="s2">&quot;OUTPUT.tif&quot;</span>,
        OpenArray<span class="o">(</span> clip, <span class="nv">prototype_ds</span><span class="o">=</span>raster_path, <span class="nv">xoff</span><span class="o">=</span>xoffset, <span class="nv">yoff</span><span class="o">=</span>yoffset <span class="o">)</span>
    <span class="o">)</span>

    <span class="c"># Save as an 8-bit jpeg for an easy, quick preview</span>
    <span class="nv">clip</span> <span class="o">=</span> clip.astype<span class="o">(</span>gdalnumeric.uint8<span class="o">)</span>
    gdalnumeric.SaveArray<span class="o">(</span>clip, <span class="s2">&quot;OUTPUT.jpg&quot;</span>, <span class="nv">format</span><span class="o">=</span><span class="s2">&quot;JPEG&quot;</span><span class="o">)</span>

    gdal.ErrorReset<span class="o">()</span>


<span class="k">if</span> <span class="nv">__name__</span> <span class="o">==</span> <span class="s1">&#39;__main__&#39;</span>:

    <span class="c">#</span>
    <span class="c"># example run : $ python clip.py /&lt;full-path&gt;/&lt;shapefile-name&gt;.shp /&lt;full-path&gt;/&lt;raster-name&gt;.tif</span>
    <span class="c">#</span>
    <span class="k">if</span> len<span class="o">(</span> sys.argv <span class="o">)</span> &lt; 2:
        print <span class="s2">&quot;[ ERROR ] you must two args. 1) the full shapefile path and 2) the full raster path&quot;</span>
        sys.exit<span class="o">(</span> <span class="m">1</span> <span class="o">)</span>

    main<span class="o">(</span> sys.argv<span class="o">[</span>1<span class="o">]</span>, sys.argv<span class="o">[</span>2<span class="o">]</span> <span class="o">)</span>
</pre></div>
</div>
<p>After Image: the clipped geotiff with the timezone border overlayed in orange on top of input geotiff:</p>
<img alt="_images/clip_raster_after.png" src="_images/clip_raster_after.png" />
</div>
<div class="section" id="calculate-zonal-statistics">
<h2>Calculate zonal statistics<a class="headerlink" href="#calculate-zonal-statistics" title="Permalink to this headline">¶</a></h2>
<p>This recipe calculates statistics on values of a raster within the zones of a vector dataset. It returns for each feature a dictionary item (FID) with the statistical values in the following order: Average, Mean, Medain, Standard Deviation, Variance</p>
<p>While this recipe works and is a good example, it is generally recommended to use [rasterstats](<a class="reference external" href="https://github.com/perrygeo/python-raster-stats">https://github.com/perrygeo/python-raster-stats</a>) for calculating zonal statistics with Python.</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">import</span> <span class="nn">gdal</span><span class="o">,</span> <span class="nn">ogr</span><span class="o">,</span> <span class="nn">osr</span><span class="o">,</span> <span class="nn">numpy</span>
<span class="kn">import</span> <span class="nn">sys</span>


<span class="k">def</span> <span class="nf">zonal_stats</span><span class="p">(</span><span class="n">feat</span><span class="p">,</span> <span class="n">input_zone_polygon</span><span class="p">,</span> <span class="n">input_value_raster</span><span class="p">):</span>

    <span class="c"># Open data</span>
    <span class="n">raster</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span><span class="n">input_value_raster</span><span class="p">)</span>
    <span class="n">shp</span> <span class="o">=</span> <span class="n">ogr</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span><span class="n">input_zone_polygon</span><span class="p">)</span>
    <span class="n">lyr</span> <span class="o">=</span> <span class="n">shp</span><span class="o">.</span><span class="n">GetLayer</span><span class="p">()</span>

    <span class="c"># Get raster georeference info</span>
    <span class="n">transform</span> <span class="o">=</span> <span class="n">raster</span><span class="o">.</span><span class="n">GetGeoTransform</span><span class="p">()</span>
    <span class="n">xOrigin</span> <span class="o">=</span> <span class="n">transform</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
    <span class="n">yOrigin</span> <span class="o">=</span> <span class="n">transform</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span>
    <span class="n">pixelWidth</span> <span class="o">=</span> <span class="n">transform</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
    <span class="n">pixelHeight</span> <span class="o">=</span> <span class="n">transform</span><span class="p">[</span><span class="mi">5</span><span class="p">]</span>

    <span class="c"># Reproject vector geometry to same projection as raster</span>
    <span class="n">sourceSR</span> <span class="o">=</span> <span class="n">lyr</span><span class="o">.</span><span class="n">GetSpatialRef</span><span class="p">()</span>
    <span class="n">targetSR</span> <span class="o">=</span> <span class="n">osr</span><span class="o">.</span><span class="n">SpatialReference</span><span class="p">()</span>
    <span class="n">targetSR</span><span class="o">.</span><span class="n">ImportFromWkt</span><span class="p">(</span><span class="n">raster</span><span class="o">.</span><span class="n">GetProjectionRef</span><span class="p">())</span>
    <span class="n">coordTrans</span> <span class="o">=</span> <span class="n">osr</span><span class="o">.</span><span class="n">CoordinateTransformation</span><span class="p">(</span><span class="n">sourceSR</span><span class="p">,</span><span class="n">targetSR</span><span class="p">)</span>
    <span class="n">feat</span> <span class="o">=</span> <span class="n">lyr</span><span class="o">.</span><span class="n">GetNextFeature</span><span class="p">()</span>
    <span class="n">geom</span> <span class="o">=</span> <span class="n">feat</span><span class="o">.</span><span class="n">GetGeometryRef</span><span class="p">()</span>
    <span class="n">geom</span><span class="o">.</span><span class="n">Transform</span><span class="p">(</span><span class="n">coordTrans</span><span class="p">)</span>

    <span class="c"># Get extent of feat</span>
    <span class="n">geom</span> <span class="o">=</span> <span class="n">feat</span><span class="o">.</span><span class="n">GetGeometryRef</span><span class="p">()</span>
    <span class="k">if</span> <span class="p">(</span><span class="n">geom</span><span class="o">.</span><span class="n">GetGeometryName</span><span class="p">()</span> <span class="o">==</span> <span class="s">&#39;MULTIPOLYGON&#39;</span><span class="p">):</span>
        <span class="n">count</span> <span class="o">=</span> <span class="mi">0</span>
        <span class="n">pointsX</span> <span class="o">=</span> <span class="p">[];</span> <span class="n">pointsY</span> <span class="o">=</span> <span class="p">[]</span>
        <span class="k">for</span> <span class="n">polygon</span> <span class="ow">in</span> <span class="n">geom</span><span class="p">:</span>
            <span class="n">geomInner</span> <span class="o">=</span> <span class="n">geom</span><span class="o">.</span><span class="n">GetGeometryRef</span><span class="p">(</span><span class="n">count</span><span class="p">)</span>
            <span class="n">ring</span> <span class="o">=</span> <span class="n">geomInner</span><span class="o">.</span><span class="n">GetGeometryRef</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span>
            <span class="n">numpoints</span> <span class="o">=</span> <span class="n">ring</span><span class="o">.</span><span class="n">GetPointCount</span><span class="p">()</span>
            <span class="k">for</span> <span class="n">p</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">numpoints</span><span class="p">):</span>
                    <span class="n">lon</span><span class="p">,</span> <span class="n">lat</span><span class="p">,</span> <span class="n">z</span> <span class="o">=</span> <span class="n">ring</span><span class="o">.</span><span class="n">GetPoint</span><span class="p">(</span><span class="n">p</span><span class="p">)</span>
                    <span class="n">pointsX</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">lon</span><span class="p">)</span>
                    <span class="n">pointsY</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">lat</span><span class="p">)</span>
            <span class="n">count</span> <span class="o">+=</span> <span class="mi">1</span>
    <span class="k">elif</span> <span class="p">(</span><span class="n">geom</span><span class="o">.</span><span class="n">GetGeometryName</span><span class="p">()</span> <span class="o">==</span> <span class="s">&#39;POLYGON&#39;</span><span class="p">):</span>
        <span class="n">ring</span> <span class="o">=</span> <span class="n">geom</span><span class="o">.</span><span class="n">GetGeometryRef</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span>
        <span class="n">numpoints</span> <span class="o">=</span> <span class="n">ring</span><span class="o">.</span><span class="n">GetPointCount</span><span class="p">()</span>
        <span class="n">pointsX</span> <span class="o">=</span> <span class="p">[];</span> <span class="n">pointsY</span> <span class="o">=</span> <span class="p">[]</span>
        <span class="k">for</span> <span class="n">p</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">numpoints</span><span class="p">):</span>
                <span class="n">lon</span><span class="p">,</span> <span class="n">lat</span><span class="p">,</span> <span class="n">z</span> <span class="o">=</span> <span class="n">ring</span><span class="o">.</span><span class="n">GetPoint</span><span class="p">(</span><span class="n">p</span><span class="p">)</span>
                <span class="n">pointsX</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">lon</span><span class="p">)</span>
                <span class="n">pointsY</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">lat</span><span class="p">)</span>

    <span class="k">else</span><span class="p">:</span>
        <span class="n">sys</span><span class="o">.</span><span class="n">exit</span><span class="p">(</span><span class="s">&quot;ERROR: Geometry needs to be either Polygon or Multipolygon&quot;</span><span class="p">)</span>

    <span class="n">xmin</span> <span class="o">=</span> <span class="nb">min</span><span class="p">(</span><span class="n">pointsX</span><span class="p">)</span>
    <span class="n">xmax</span> <span class="o">=</span> <span class="nb">max</span><span class="p">(</span><span class="n">pointsX</span><span class="p">)</span>
    <span class="n">ymin</span> <span class="o">=</span> <span class="nb">min</span><span class="p">(</span><span class="n">pointsY</span><span class="p">)</span>
    <span class="n">ymax</span> <span class="o">=</span> <span class="nb">max</span><span class="p">(</span><span class="n">pointsY</span><span class="p">)</span>

    <span class="c"># Specify offset and rows and columns to read</span>
    <span class="n">xoff</span> <span class="o">=</span> <span class="nb">int</span><span class="p">((</span><span class="n">xmin</span> <span class="o">-</span> <span class="n">xOrigin</span><span class="p">)</span><span class="o">/</span><span class="n">pixelWidth</span><span class="p">)</span>
    <span class="n">yoff</span> <span class="o">=</span> <span class="nb">int</span><span class="p">((</span><span class="n">yOrigin</span> <span class="o">-</span> <span class="n">ymax</span><span class="p">)</span><span class="o">/</span><span class="n">pixelWidth</span><span class="p">)</span>
    <span class="n">xcount</span> <span class="o">=</span> <span class="nb">int</span><span class="p">((</span><span class="n">xmax</span> <span class="o">-</span> <span class="n">xmin</span><span class="p">)</span><span class="o">/</span><span class="n">pixelWidth</span><span class="p">)</span><span class="o">+</span><span class="mi">1</span>
    <span class="n">ycount</span> <span class="o">=</span> <span class="nb">int</span><span class="p">((</span><span class="n">ymax</span> <span class="o">-</span> <span class="n">ymin</span><span class="p">)</span><span class="o">/</span><span class="n">pixelWidth</span><span class="p">)</span><span class="o">+</span><span class="mi">1</span>

    <span class="c"># Create memory target raster</span>
    <span class="n">target_ds</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">GetDriverByName</span><span class="p">(</span><span class="s">&#39;MEM&#39;</span><span class="p">)</span><span class="o">.</span><span class="n">Create</span><span class="p">(</span><span class="s">&#39;&#39;</span><span class="p">,</span> <span class="n">xcount</span><span class="p">,</span> <span class="n">ycount</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="n">gdal</span><span class="o">.</span><span class="n">GDT_Byte</span><span class="p">)</span>
    <span class="n">target_ds</span><span class="o">.</span><span class="n">SetGeoTransform</span><span class="p">((</span>
        <span class="n">xmin</span><span class="p">,</span> <span class="n">pixelWidth</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span>
        <span class="n">ymax</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="n">pixelHeight</span><span class="p">,</span>
    <span class="p">))</span>

    <span class="c"># Create for target raster the same projection as for the value raster</span>
    <span class="n">raster_srs</span> <span class="o">=</span> <span class="n">osr</span><span class="o">.</span><span class="n">SpatialReference</span><span class="p">()</span>
    <span class="n">raster_srs</span><span class="o">.</span><span class="n">ImportFromWkt</span><span class="p">(</span><span class="n">raster</span><span class="o">.</span><span class="n">GetProjectionRef</span><span class="p">())</span>
    <span class="n">target_ds</span><span class="o">.</span><span class="n">SetProjection</span><span class="p">(</span><span class="n">raster_srs</span><span class="o">.</span><span class="n">ExportToWkt</span><span class="p">())</span>

    <span class="c"># Rasterize zone polygon to raster</span>
    <span class="n">gdal</span><span class="o">.</span><span class="n">RasterizeLayer</span><span class="p">(</span><span class="n">target_ds</span><span class="p">,</span> <span class="p">[</span><span class="mi">1</span><span class="p">],</span> <span class="n">lyr</span><span class="p">,</span> <span class="n">burn_values</span><span class="o">=</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>

    <span class="c"># Read raster as arrays</span>
    <span class="n">banddataraster</span> <span class="o">=</span> <span class="n">raster</span><span class="o">.</span><span class="n">GetRasterBand</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>
    <span class="n">dataraster</span> <span class="o">=</span> <span class="n">banddataraster</span><span class="o">.</span><span class="n">ReadAsArray</span><span class="p">(</span><span class="n">xoff</span><span class="p">,</span> <span class="n">yoff</span><span class="p">,</span> <span class="n">xcount</span><span class="p">,</span> <span class="n">ycount</span><span class="p">)</span><span class="o">.</span><span class="n">astype</span><span class="p">(</span><span class="n">numpy</span><span class="o">.</span><span class="n">float</span><span class="p">)</span>

    <span class="n">bandmask</span> <span class="o">=</span> <span class="n">target_ds</span><span class="o">.</span><span class="n">GetRasterBand</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>
    <span class="n">datamask</span> <span class="o">=</span> <span class="n">bandmask</span><span class="o">.</span><span class="n">ReadAsArray</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="n">xcount</span><span class="p">,</span> <span class="n">ycount</span><span class="p">)</span><span class="o">.</span><span class="n">astype</span><span class="p">(</span><span class="n">numpy</span><span class="o">.</span><span class="n">float</span><span class="p">)</span>

    <span class="c"># Mask zone of raster</span>
    <span class="n">zoneraster</span> <span class="o">=</span> <span class="n">numpy</span><span class="o">.</span><span class="n">ma</span><span class="o">.</span><span class="n">masked_array</span><span class="p">(</span><span class="n">dataraster</span><span class="p">,</span>  <span class="n">numpy</span><span class="o">.</span><span class="n">logical_not</span><span class="p">(</span><span class="n">datamask</span><span class="p">))</span>

    <span class="c"># Calculate statistics of zonal raster</span>
    <span class="k">return</span> <span class="n">numpy</span><span class="o">.</span><span class="n">average</span><span class="p">(</span><span class="n">zoneraster</span><span class="p">),</span><span class="n">numpy</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">zoneraster</span><span class="p">),</span><span class="n">numpy</span><span class="o">.</span><span class="n">median</span><span class="p">(</span><span class="n">zoneraster</span><span class="p">),</span><span class="n">numpy</span><span class="o">.</span><span class="n">std</span><span class="p">(</span><span class="n">zoneraster</span><span class="p">),</span><span class="n">numpy</span><span class="o">.</span><span class="n">var</span><span class="p">(</span><span class="n">zoneraster</span><span class="p">)</span>


<span class="k">def</span> <span class="nf">loop_zonal_stats</span><span class="p">(</span><span class="n">input_zone_polygon</span><span class="p">,</span> <span class="n">input_value_raster</span><span class="p">):</span>

    <span class="n">shp</span> <span class="o">=</span> <span class="n">ogr</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span><span class="n">input_zone_polygon</span><span class="p">)</span>
    <span class="n">lyr</span> <span class="o">=</span> <span class="n">shp</span><span class="o">.</span><span class="n">GetLayer</span><span class="p">()</span>
    <span class="n">featList</span> <span class="o">=</span> <span class="nb">range</span><span class="p">(</span><span class="n">lyr</span><span class="o">.</span><span class="n">GetFeatureCount</span><span class="p">())</span>
    <span class="n">statDict</span> <span class="o">=</span> <span class="p">{}</span>

    <span class="k">for</span> <span class="n">FID</span> <span class="ow">in</span> <span class="n">featList</span><span class="p">:</span>
        <span class="n">feat</span> <span class="o">=</span> <span class="n">lyr</span><span class="o">.</span><span class="n">GetFeature</span><span class="p">(</span><span class="n">FID</span><span class="p">)</span>
        <span class="n">meanValue</span> <span class="o">=</span> <span class="n">zonal_stats</span><span class="p">(</span><span class="n">feat</span><span class="p">,</span> <span class="n">input_zone_polygon</span><span class="p">,</span> <span class="n">input_value_raster</span><span class="p">)</span>
        <span class="n">statDict</span><span class="p">[</span><span class="n">FID</span><span class="p">]</span> <span class="o">=</span> <span class="n">meanValue</span>
    <span class="k">return</span> <span class="n">statDict</span>

<span class="k">def</span> <span class="nf">main</span><span class="p">(</span><span class="n">input_zone_polygon</span><span class="p">,</span> <span class="n">input_value_raster</span><span class="p">):</span>
    <span class="k">return</span> <span class="n">loop_zonal_stats</span><span class="p">(</span><span class="n">input_zone_polygon</span><span class="p">,</span> <span class="n">input_value_raster</span><span class="p">)</span>


<span class="k">if</span> <span class="n">__name__</span> <span class="o">==</span> <span class="s">&quot;__main__&quot;</span><span class="p">:</span>

    <span class="c">#</span>
    <span class="c"># Returns for each feature a dictionary item (FID) with the statistical values in the following order: Average, Mean, Medain, Standard Deviation, Variance</span>
    <span class="c">#</span>
    <span class="c"># example run : $ python grid.py &lt;full-path&gt;&lt;output-shapefile-name&gt;.shp xmin xmax ymin ymax gridHeight gridWidth</span>
    <span class="c">#</span>

    <span class="k">if</span> <span class="nb">len</span><span class="p">(</span> <span class="n">sys</span><span class="o">.</span><span class="n">argv</span> <span class="p">)</span> <span class="o">!=</span> <span class="mi">3</span><span class="p">:</span>
        <span class="k">print</span> <span class="s">&quot;[ ERROR ] you must supply two arguments: input-zone-shapefile-name.shp input-value-raster-name.tif &quot;</span>
        <span class="n">sys</span><span class="o">.</span><span class="n">exit</span><span class="p">(</span> <span class="mi">1</span> <span class="p">)</span>
    <span class="k">print</span> <span class="s">&#39;Returns for each feature a dictionary item (FID) with the statistical values in the following order: Average, Mean, Medain, Standard Deviation, Variance&#39;</span>
    <span class="k">print</span> <span class="n">main</span><span class="p">(</span> <span class="n">sys</span><span class="o">.</span><span class="n">argv</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span> <span class="n">sys</span><span class="o">.</span><span class="n">argv</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span> <span class="p">)</span>
</pre></div>
</div>
</div>
<div class="section" id="raster-to-vector-line">
<h2>Raster to vector line<a class="headerlink" href="#raster-to-vector-line" title="Permalink to this headline">¶</a></h2>
<p>This recipe converts raster pixels with a specified value to vector lines. For example the blue pixels (value = 0) are converted to vector lines.</p>
<img alt="_images/raster2line1.png" src="_images/raster2line1.png" />
<div class="highlight-python"><div class="highlight"><pre><span class="kn">import</span> <span class="nn">ogr</span><span class="o">,</span> <span class="nn">gdal</span><span class="o">,</span> <span class="nn">osr</span><span class="o">,</span> <span class="nn">os</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="kn">as</span> <span class="nn">np</span>
<span class="kn">import</span> <span class="nn">itertools</span>
<span class="kn">from</span> <span class="nn">math</span> <span class="kn">import</span> <span class="n">sqrt</span><span class="p">,</span><span class="n">ceil</span>

<span class="k">def</span> <span class="nf">pixelOffset2coord</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">,</span><span class="n">xOffset</span><span class="p">,</span><span class="n">yOffset</span><span class="p">):</span>
    <span class="n">raster</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">)</span>
    <span class="n">geotransform</span> <span class="o">=</span> <span class="n">raster</span><span class="o">.</span><span class="n">GetGeoTransform</span><span class="p">()</span>
    <span class="n">originX</span> <span class="o">=</span> <span class="n">geotransform</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
    <span class="n">originY</span> <span class="o">=</span> <span class="n">geotransform</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span>
    <span class="n">pixelWidth</span> <span class="o">=</span> <span class="n">geotransform</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
    <span class="n">pixelHeight</span> <span class="o">=</span> <span class="n">geotransform</span><span class="p">[</span><span class="mi">5</span><span class="p">]</span>
    <span class="n">coordX</span> <span class="o">=</span> <span class="n">originX</span><span class="o">+</span><span class="n">pixelWidth</span><span class="o">*</span><span class="n">xOffset</span>
    <span class="n">coordY</span> <span class="o">=</span> <span class="n">originY</span><span class="o">+</span><span class="n">pixelHeight</span><span class="o">*</span><span class="n">yOffset</span>
    <span class="k">return</span> <span class="n">coordX</span><span class="p">,</span> <span class="n">coordY</span>

<span class="k">def</span> <span class="nf">raster2array</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">):</span>
    <span class="n">raster</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">)</span>
    <span class="n">band</span> <span class="o">=</span> <span class="n">raster</span><span class="o">.</span><span class="n">GetRasterBand</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>
    <span class="n">array</span> <span class="o">=</span> <span class="n">band</span><span class="o">.</span><span class="n">ReadAsArray</span><span class="p">()</span>
    <span class="k">return</span> <span class="n">array</span>

<span class="k">def</span> <span class="nf">array2shp</span><span class="p">(</span><span class="n">array</span><span class="p">,</span><span class="n">outSHPfn</span><span class="p">,</span><span class="n">rasterfn</span><span class="p">,</span><span class="n">pixelValue</span><span class="p">):</span>

    <span class="c"># max distance between points</span>
    <span class="n">raster</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">)</span>
    <span class="n">geotransform</span> <span class="o">=</span> <span class="n">raster</span><span class="o">.</span><span class="n">GetGeoTransform</span><span class="p">()</span>
    <span class="n">pixelWidth</span> <span class="o">=</span> <span class="n">geotransform</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
    <span class="n">maxDistance</span> <span class="o">=</span> <span class="n">ceil</span><span class="p">(</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">pixelWidth</span><span class="o">*</span><span class="n">pixelWidth</span><span class="p">))</span>
    <span class="k">print</span> <span class="n">maxDistance</span>

    <span class="c"># array2dict</span>
    <span class="n">count</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="n">roadList</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">where</span><span class="p">(</span><span class="n">array</span> <span class="o">==</span> <span class="n">pixelValue</span><span class="p">)</span>
    <span class="n">multipoint</span> <span class="o">=</span> <span class="n">ogr</span><span class="o">.</span><span class="n">Geometry</span><span class="p">(</span><span class="n">ogr</span><span class="o">.</span><span class="n">wkbMultiLineString</span><span class="p">)</span>
    <span class="n">pointDict</span> <span class="o">=</span> <span class="p">{}</span>
    <span class="k">for</span> <span class="n">indexY</span> <span class="ow">in</span> <span class="n">roadList</span><span class="p">[</span><span class="mi">0</span><span class="p">]:</span>
        <span class="n">indexX</span> <span class="o">=</span> <span class="n">roadList</span><span class="p">[</span><span class="mi">1</span><span class="p">][</span><span class="n">count</span><span class="p">]</span>
        <span class="n">Xcoord</span><span class="p">,</span> <span class="n">Ycoord</span> <span class="o">=</span> <span class="n">pixelOffset2coord</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">,</span><span class="n">indexX</span><span class="p">,</span><span class="n">indexY</span><span class="p">)</span>
        <span class="n">pointDict</span><span class="p">[</span><span class="n">count</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">Xcoord</span><span class="p">,</span> <span class="n">Ycoord</span><span class="p">)</span>
        <span class="n">count</span> <span class="o">+=</span> <span class="mi">1</span>

    <span class="c"># dict2wkbMultiLineString</span>
    <span class="n">multiline</span> <span class="o">=</span> <span class="n">ogr</span><span class="o">.</span><span class="n">Geometry</span><span class="p">(</span><span class="n">ogr</span><span class="o">.</span><span class="n">wkbMultiLineString</span><span class="p">)</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">itertools</span><span class="o">.</span><span class="n">combinations</span><span class="p">(</span><span class="n">pointDict</span><span class="o">.</span><span class="n">values</span><span class="p">(),</span> <span class="mi">2</span><span class="p">):</span>
        <span class="n">point1</span> <span class="o">=</span> <span class="n">ogr</span><span class="o">.</span><span class="n">Geometry</span><span class="p">(</span><span class="n">ogr</span><span class="o">.</span><span class="n">wkbPoint</span><span class="p">)</span>
        <span class="n">point1</span><span class="o">.</span><span class="n">AddPoint</span><span class="p">(</span><span class="n">i</span><span class="p">[</span><span class="mi">0</span><span class="p">][</span><span class="mi">0</span><span class="p">],</span><span class="n">i</span><span class="p">[</span><span class="mi">0</span><span class="p">][</span><span class="mi">1</span><span class="p">])</span>
        <span class="n">point2</span> <span class="o">=</span> <span class="n">ogr</span><span class="o">.</span><span class="n">Geometry</span><span class="p">(</span><span class="n">ogr</span><span class="o">.</span><span class="n">wkbPoint</span><span class="p">)</span>
        <span class="n">point2</span><span class="o">.</span><span class="n">AddPoint</span><span class="p">(</span><span class="n">i</span><span class="p">[</span><span class="mi">1</span><span class="p">][</span><span class="mi">0</span><span class="p">],</span><span class="n">i</span><span class="p">[</span><span class="mi">1</span><span class="p">][</span><span class="mi">1</span><span class="p">])</span>

        <span class="n">distance</span> <span class="o">=</span> <span class="n">point1</span><span class="o">.</span><span class="n">Distance</span><span class="p">(</span><span class="n">point2</span><span class="p">)</span>

        <span class="k">if</span> <span class="n">distance</span> <span class="o">&lt;</span> <span class="n">maxDistance</span><span class="p">:</span>
            <span class="n">line</span> <span class="o">=</span> <span class="n">ogr</span><span class="o">.</span><span class="n">Geometry</span><span class="p">(</span><span class="n">ogr</span><span class="o">.</span><span class="n">wkbLineString</span><span class="p">)</span>
            <span class="n">line</span><span class="o">.</span><span class="n">AddPoint</span><span class="p">(</span><span class="n">i</span><span class="p">[</span><span class="mi">0</span><span class="p">][</span><span class="mi">0</span><span class="p">],</span><span class="n">i</span><span class="p">[</span><span class="mi">0</span><span class="p">][</span><span class="mi">1</span><span class="p">])</span>
            <span class="n">line</span><span class="o">.</span><span class="n">AddPoint</span><span class="p">(</span><span class="n">i</span><span class="p">[</span><span class="mi">1</span><span class="p">][</span><span class="mi">0</span><span class="p">],</span><span class="n">i</span><span class="p">[</span><span class="mi">1</span><span class="p">][</span><span class="mi">1</span><span class="p">])</span>
            <span class="n">multiline</span><span class="o">.</span><span class="n">AddGeometry</span><span class="p">(</span><span class="n">line</span><span class="p">)</span>

    <span class="c"># wkbMultiLineString2shp</span>
    <span class="n">shpDriver</span> <span class="o">=</span> <span class="n">ogr</span><span class="o">.</span><span class="n">GetDriverByName</span><span class="p">(</span><span class="s">&quot;ESRI Shapefile&quot;</span><span class="p">)</span>
    <span class="k">if</span> <span class="n">os</span><span class="o">.</span><span class="n">path</span><span class="o">.</span><span class="n">exists</span><span class="p">(</span><span class="n">outSHPfn</span><span class="p">):</span>
        <span class="n">shpDriver</span><span class="o">.</span><span class="n">DeleteDataSource</span><span class="p">(</span><span class="n">outSHPfn</span><span class="p">)</span>
    <span class="n">outDataSource</span> <span class="o">=</span> <span class="n">shpDriver</span><span class="o">.</span><span class="n">CreateDataSource</span><span class="p">(</span><span class="n">outSHPfn</span><span class="p">)</span>
    <span class="n">outLayer</span> <span class="o">=</span> <span class="n">outDataSource</span><span class="o">.</span><span class="n">CreateLayer</span><span class="p">(</span><span class="n">outSHPfn</span><span class="p">,</span> <span class="n">geom_type</span><span class="o">=</span><span class="n">ogr</span><span class="o">.</span><span class="n">wkbMultiLineString</span> <span class="p">)</span>
    <span class="n">featureDefn</span> <span class="o">=</span> <span class="n">outLayer</span><span class="o">.</span><span class="n">GetLayerDefn</span><span class="p">()</span>
    <span class="n">outFeature</span> <span class="o">=</span> <span class="n">ogr</span><span class="o">.</span><span class="n">Feature</span><span class="p">(</span><span class="n">featureDefn</span><span class="p">)</span>
    <span class="n">outFeature</span><span class="o">.</span><span class="n">SetGeometry</span><span class="p">(</span><span class="n">multiline</span><span class="p">)</span>
    <span class="n">outLayer</span><span class="o">.</span><span class="n">CreateFeature</span><span class="p">(</span><span class="n">outFeature</span><span class="p">)</span>


<span class="k">def</span> <span class="nf">main</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">,</span><span class="n">outSHPfn</span><span class="p">,</span><span class="n">pixelValue</span><span class="p">):</span>
    <span class="n">array</span> <span class="o">=</span> <span class="n">raster2array</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">)</span>
    <span class="n">array2shp</span><span class="p">(</span><span class="n">array</span><span class="p">,</span><span class="n">outSHPfn</span><span class="p">,</span><span class="n">rasterfn</span><span class="p">,</span><span class="n">pixelValue</span><span class="p">)</span>

<span class="k">if</span> <span class="n">__name__</span> <span class="o">==</span> <span class="s">&quot;__main__&quot;</span><span class="p">:</span>
    <span class="n">rasterfn</span> <span class="o">=</span> <span class="s">&#39;test.tif&#39;</span>
    <span class="n">outSHPfn</span> <span class="o">=</span> <span class="s">&#39;test.shp&#39;</span>
    <span class="n">pixelValue</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="n">main</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">,</span><span class="n">outSHPfn</span><span class="p">,</span><span class="n">pixelValue</span><span class="p">)</span>
</pre></div>
</div>
<p>The green lines show the converted lines.</p>
<img alt="_images/raster2line2.png" src="_images/raster2line2.png" />
</div>
<div class="section" id="create-raster-from-array">
<h2>Create raster from array<a class="headerlink" href="#create-raster-from-array" title="Permalink to this headline">¶</a></h2>
<p>This recipe creates a raster from an array</p>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">import</span> <span class="nn">gdal</span><span class="o">,</span> <span class="nn">ogr</span><span class="o">,</span> <span class="nn">os</span><span class="o">,</span> <span class="nn">osr</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="kn">as</span> <span class="nn">np</span>


<span class="k">def</span> <span class="nf">array2raster</span><span class="p">(</span><span class="n">newRasterfn</span><span class="p">,</span><span class="n">rasterOrigin</span><span class="p">,</span><span class="n">pixelWidth</span><span class="p">,</span><span class="n">pixelHeight</span><span class="p">,</span><span class="n">array</span><span class="p">):</span>

    <span class="n">cols</span> <span class="o">=</span> <span class="n">array</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
    <span class="n">rows</span> <span class="o">=</span> <span class="n">array</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
    <span class="n">originX</span> <span class="o">=</span> <span class="n">rasterOrigin</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
    <span class="n">originY</span> <span class="o">=</span> <span class="n">rasterOrigin</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>

    <span class="n">driver</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">GetDriverByName</span><span class="p">(</span><span class="s">&#39;GTiff&#39;</span><span class="p">)</span>
    <span class="n">outRaster</span> <span class="o">=</span> <span class="n">driver</span><span class="o">.</span><span class="n">Create</span><span class="p">(</span><span class="n">newRasterfn</span><span class="p">,</span> <span class="n">cols</span><span class="p">,</span> <span class="n">rows</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="n">gdal</span><span class="o">.</span><span class="n">GDT_Byte</span><span class="p">)</span>
    <span class="n">outRaster</span><span class="o">.</span><span class="n">SetGeoTransform</span><span class="p">((</span><span class="n">originX</span><span class="p">,</span> <span class="n">pixelWidth</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="n">originY</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="n">pixelHeight</span><span class="p">))</span>
    <span class="n">outband</span> <span class="o">=</span> <span class="n">outRaster</span><span class="o">.</span><span class="n">GetRasterBand</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>
    <span class="n">outband</span><span class="o">.</span><span class="n">WriteArray</span><span class="p">(</span><span class="n">array</span><span class="p">)</span>
    <span class="n">outRasterSRS</span> <span class="o">=</span> <span class="n">osr</span><span class="o">.</span><span class="n">SpatialReference</span><span class="p">()</span>
    <span class="n">outRasterSRS</span><span class="o">.</span><span class="n">ImportFromEPSG</span><span class="p">(</span><span class="mi">4326</span><span class="p">)</span>
    <span class="n">outRaster</span><span class="o">.</span><span class="n">SetProjection</span><span class="p">(</span><span class="n">outRasterSRS</span><span class="o">.</span><span class="n">ExportToWkt</span><span class="p">())</span>
    <span class="n">outband</span><span class="o">.</span><span class="n">FlushCache</span><span class="p">()</span>


<span class="k">def</span> <span class="nf">main</span><span class="p">(</span><span class="n">newRasterfn</span><span class="p">,</span><span class="n">rasterOrigin</span><span class="p">,</span><span class="n">pixelWidth</span><span class="p">,</span><span class="n">pixelHeight</span><span class="p">,</span><span class="n">array</span><span class="p">):</span>
    <span class="n">reversed_arr</span> <span class="o">=</span> <span class="n">array</span><span class="p">[::</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="c"># reverse array so the tif looks like the array</span>
    <span class="n">array2raster</span><span class="p">(</span><span class="n">newRasterfn</span><span class="p">,</span><span class="n">rasterOrigin</span><span class="p">,</span><span class="n">pixelWidth</span><span class="p">,</span><span class="n">pixelHeight</span><span class="p">,</span><span class="n">reversed_arr</span><span class="p">)</span> <span class="c"># convert array to raster</span>


<span class="k">if</span> <span class="n">__name__</span> <span class="o">==</span> <span class="s">&quot;__main__&quot;</span><span class="p">:</span>
    <span class="n">rasterOrigin</span> <span class="o">=</span> <span class="p">(</span><span class="o">-</span><span class="mf">123.25745</span><span class="p">,</span><span class="mf">45.43013</span><span class="p">)</span>
    <span class="n">pixelWidth</span> <span class="o">=</span> <span class="mi">10</span>
    <span class="n">pixelHeight</span> <span class="o">=</span> <span class="mi">10</span>
    <span class="n">newRasterfn</span> <span class="o">=</span> <span class="s">&#39;test.tif&#39;</span>
    <span class="n">array</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([[</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span>
                      <span class="p">[</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span>
                      <span class="p">[</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span>
                      <span class="p">[</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span>
                      <span class="p">[</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span>
                      <span class="p">[</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span>
                      <span class="p">[</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span>
                      <span class="p">[</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span>
                      <span class="p">[</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span>
                      <span class="p">[</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">]])</span>


    <span class="n">main</span><span class="p">(</span><span class="n">newRasterfn</span><span class="p">,</span><span class="n">rasterOrigin</span><span class="p">,</span><span class="n">pixelWidth</span><span class="p">,</span><span class="n">pixelHeight</span><span class="p">,</span><span class="n">array</span><span class="p">)</span>
</pre></div>
</div>
<p>The create raster from the example array looks like this:</p>
<img alt="_images/array2raster.png" src="_images/array2raster.png" />
</div>
<div class="section" id="create-least-cost-path">
<h2>Create least cost path<a class="headerlink" href="#create-least-cost-path" title="Permalink to this headline">¶</a></h2>
<p>This recipe creates a least cost path between two coordinates based on a raster cost surface.
In the example below, a cost path between point 1 and point 2 is created based on a slope raster.</p>
<img alt="_images/costpath1.png" src="_images/costpath1.png" />
<div class="highlight-python"><div class="highlight"><pre><span class="kn">import</span> <span class="nn">gdal</span><span class="o">,</span> <span class="nn">osr</span>
<span class="kn">from</span> <span class="nn">skimage.graph</span> <span class="kn">import</span> <span class="n">route_through_array</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="kn">as</span> <span class="nn">np</span>


<span class="k">def</span> <span class="nf">raster2array</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">):</span>
    <span class="n">raster</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">)</span>
    <span class="n">band</span> <span class="o">=</span> <span class="n">raster</span><span class="o">.</span><span class="n">GetRasterBand</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>
    <span class="n">array</span> <span class="o">=</span> <span class="n">band</span><span class="o">.</span><span class="n">ReadAsArray</span><span class="p">()</span>
    <span class="k">return</span> <span class="n">array</span>

<span class="k">def</span> <span class="nf">coord2pixelOffset</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">,</span><span class="n">x</span><span class="p">,</span><span class="n">y</span><span class="p">):</span>
    <span class="n">raster</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">)</span>
    <span class="n">geotransform</span> <span class="o">=</span> <span class="n">raster</span><span class="o">.</span><span class="n">GetGeoTransform</span><span class="p">()</span>
    <span class="n">originX</span> <span class="o">=</span> <span class="n">geotransform</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
    <span class="n">originY</span> <span class="o">=</span> <span class="n">geotransform</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span>
    <span class="n">pixelWidth</span> <span class="o">=</span> <span class="n">geotransform</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
    <span class="n">pixelHeight</span> <span class="o">=</span> <span class="n">geotransform</span><span class="p">[</span><span class="mi">5</span><span class="p">]</span>
    <span class="n">xOffset</span> <span class="o">=</span> <span class="nb">int</span><span class="p">((</span><span class="n">x</span> <span class="o">-</span> <span class="n">originX</span><span class="p">)</span><span class="o">/</span><span class="n">pixelWidth</span><span class="p">)</span>
    <span class="n">yOffset</span> <span class="o">=</span> <span class="nb">int</span><span class="p">((</span><span class="n">y</span> <span class="o">-</span> <span class="n">originY</span><span class="p">)</span><span class="o">/</span><span class="n">pixelHeight</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">xOffset</span><span class="p">,</span><span class="n">yOffset</span>

<span class="k">def</span> <span class="nf">createPath</span><span class="p">(</span><span class="n">CostSurfacefn</span><span class="p">,</span><span class="n">costSurfaceArray</span><span class="p">,</span><span class="n">startCoord</span><span class="p">,</span><span class="n">stopCoord</span><span class="p">):</span>

    <span class="c"># coordinates to array index</span>
    <span class="n">startCoordX</span> <span class="o">=</span> <span class="n">startCoord</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
    <span class="n">startCoordY</span> <span class="o">=</span> <span class="n">startCoord</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
    <span class="n">startIndexX</span><span class="p">,</span><span class="n">startIndexY</span> <span class="o">=</span> <span class="n">coord2pixelOffset</span><span class="p">(</span><span class="n">CostSurfacefn</span><span class="p">,</span><span class="n">startCoordX</span><span class="p">,</span><span class="n">startCoordY</span><span class="p">)</span>

    <span class="n">stopCoordX</span> <span class="o">=</span> <span class="n">stopCoord</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
    <span class="n">stopCoordY</span> <span class="o">=</span> <span class="n">stopCoord</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
    <span class="n">stopIndexX</span><span class="p">,</span><span class="n">stopIndexY</span> <span class="o">=</span> <span class="n">coord2pixelOffset</span><span class="p">(</span><span class="n">CostSurfacefn</span><span class="p">,</span><span class="n">stopCoordX</span><span class="p">,</span><span class="n">stopCoordY</span><span class="p">)</span>

    <span class="c"># create path</span>
    <span class="n">indices</span><span class="p">,</span> <span class="n">weight</span> <span class="o">=</span> <span class="n">route_through_array</span><span class="p">(</span><span class="n">costSurfaceArray</span><span class="p">,</span> <span class="p">(</span><span class="n">startIndexY</span><span class="p">,</span><span class="n">startIndexX</span><span class="p">),</span> <span class="p">(</span><span class="n">stopIndexY</span><span class="p">,</span><span class="n">stopIndexX</span><span class="p">),</span><span class="n">geometric</span><span class="o">=</span><span class="bp">True</span><span class="p">,</span><span class="n">fully_connected</span><span class="o">=</span><span class="bp">True</span><span class="p">)</span>
    <span class="n">indices</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">indices</span><span class="p">)</span><span class="o">.</span><span class="n">T</span>
    <span class="n">path</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros_like</span><span class="p">(</span><span class="n">costSurfaceArray</span><span class="p">)</span>
    <span class="n">path</span><span class="p">[</span><span class="n">indices</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span> <span class="n">indices</span><span class="p">[</span><span class="mi">1</span><span class="p">]]</span> <span class="o">=</span> <span class="mi">1</span>
    <span class="k">return</span> <span class="n">path</span>

<span class="k">def</span> <span class="nf">array2raster</span><span class="p">(</span><span class="n">newRasterfn</span><span class="p">,</span><span class="n">rasterfn</span><span class="p">,</span><span class="n">array</span><span class="p">):</span>
    <span class="n">raster</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">)</span>
    <span class="n">geotransform</span> <span class="o">=</span> <span class="n">raster</span><span class="o">.</span><span class="n">GetGeoTransform</span><span class="p">()</span>
    <span class="n">originX</span> <span class="o">=</span> <span class="n">geotransform</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
    <span class="n">originY</span> <span class="o">=</span> <span class="n">geotransform</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span>
    <span class="n">pixelWidth</span> <span class="o">=</span> <span class="n">geotransform</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
    <span class="n">pixelHeight</span> <span class="o">=</span> <span class="n">geotransform</span><span class="p">[</span><span class="mi">5</span><span class="p">]</span>
    <span class="n">cols</span> <span class="o">=</span> <span class="n">array</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
    <span class="n">rows</span> <span class="o">=</span> <span class="n">array</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>

    <span class="n">driver</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">GetDriverByName</span><span class="p">(</span><span class="s">&#39;GTiff&#39;</span><span class="p">)</span>
    <span class="n">outRaster</span> <span class="o">=</span> <span class="n">driver</span><span class="o">.</span><span class="n">Create</span><span class="p">(</span><span class="n">newRasterfn</span><span class="p">,</span> <span class="n">cols</span><span class="p">,</span> <span class="n">rows</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="n">gdal</span><span class="o">.</span><span class="n">GDT_Byte</span><span class="p">)</span>
    <span class="n">outRaster</span><span class="o">.</span><span class="n">SetGeoTransform</span><span class="p">((</span><span class="n">originX</span><span class="p">,</span> <span class="n">pixelWidth</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="n">originY</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="n">pixelHeight</span><span class="p">))</span>
    <span class="n">outband</span> <span class="o">=</span> <span class="n">outRaster</span><span class="o">.</span><span class="n">GetRasterBand</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>
    <span class="n">outband</span><span class="o">.</span><span class="n">WriteArray</span><span class="p">(</span><span class="n">array</span><span class="p">)</span>
    <span class="n">outRasterSRS</span> <span class="o">=</span> <span class="n">osr</span><span class="o">.</span><span class="n">SpatialReference</span><span class="p">()</span>
    <span class="n">outRasterSRS</span><span class="o">.</span><span class="n">ImportFromWkt</span><span class="p">(</span><span class="n">raster</span><span class="o">.</span><span class="n">GetProjectionRef</span><span class="p">())</span>
    <span class="n">outRaster</span><span class="o">.</span><span class="n">SetProjection</span><span class="p">(</span><span class="n">outRasterSRS</span><span class="o">.</span><span class="n">ExportToWkt</span><span class="p">())</span>
    <span class="n">outband</span><span class="o">.</span><span class="n">FlushCache</span><span class="p">()</span>

<span class="k">def</span> <span class="nf">main</span><span class="p">(</span><span class="n">CostSurfacefn</span><span class="p">,</span><span class="n">outputPathfn</span><span class="p">,</span><span class="n">startCoord</span><span class="p">,</span><span class="n">stopCoord</span><span class="p">):</span>

    <span class="n">costSurfaceArray</span> <span class="o">=</span> <span class="n">raster2array</span><span class="p">(</span><span class="n">CostSurfacefn</span><span class="p">)</span> <span class="c"># creates array from cost surface raster</span>

    <span class="n">pathArray</span> <span class="o">=</span> <span class="n">createPath</span><span class="p">(</span><span class="n">CostSurfacefn</span><span class="p">,</span><span class="n">costSurfaceArray</span><span class="p">,</span><span class="n">startCoord</span><span class="p">,</span><span class="n">stopCoord</span><span class="p">)</span> <span class="c"># creates path array</span>

    <span class="n">array2raster</span><span class="p">(</span><span class="n">outputPathfn</span><span class="p">,</span><span class="n">CostSurfacefn</span><span class="p">,</span><span class="n">pathArray</span><span class="p">)</span> <span class="c"># converts path array to raster</span>


<span class="k">if</span> <span class="n">__name__</span> <span class="o">==</span> <span class="s">&quot;__main__&quot;</span><span class="p">:</span>
    <span class="n">CostSurfacefn</span> <span class="o">=</span> <span class="s">&#39;CostSurface.tif&#39;</span>
    <span class="n">startCoord</span> <span class="o">=</span> <span class="p">(</span><span class="mf">345387.871</span><span class="p">,</span><span class="mf">1267855.277</span><span class="p">)</span>
    <span class="n">stopCoord</span> <span class="o">=</span> <span class="p">(</span><span class="mf">345479.425</span><span class="p">,</span><span class="mf">1267799.626</span><span class="p">)</span>
    <span class="n">outputPathfn</span> <span class="o">=</span> <span class="s">&#39;Path.tif&#39;</span>
    <span class="n">main</span><span class="p">(</span><span class="n">CostSurfacefn</span><span class="p">,</span><span class="n">outputPathfn</span><span class="p">,</span><span class="n">startCoord</span><span class="p">,</span><span class="n">stopCoord</span><span class="p">)</span>
</pre></div>
</div>
<p>The created cost path between the points.</p>
<img alt="_images/costpath2.png" src="_images/costpath2.png" />
</div>
<div class="section" id="replace-no-data-value-of-raster-with-new-value">
<h2>Replace No Data Value of Raster with new value<a class="headerlink" href="#replace-no-data-value-of-raster-with-new-value" title="Permalink to this headline">¶</a></h2>
<p>This recipe replaces the No Data Value with a new specified value</p>
<img alt="_images/costpath1.png" src="_images/costpath1.png" />
<div class="highlight-python"><div class="highlight"><pre><span class="kn">import</span> <span class="nn">gdal</span><span class="o">,</span> <span class="nn">ogr</span><span class="o">,</span> <span class="nn">osr</span><span class="o">,</span> <span class="nn">os</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="kn">as</span> <span class="nn">np</span>

<span class="k">def</span> <span class="nf">raster2array</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">):</span>
    <span class="n">raster</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">)</span>
    <span class="n">band</span> <span class="o">=</span> <span class="n">raster</span><span class="o">.</span><span class="n">GetRasterBand</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">band</span><span class="o">.</span><span class="n">ReadAsArray</span><span class="p">()</span>

<span class="k">def</span> <span class="nf">getNoDataValue</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">):</span>
    <span class="n">raster</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">)</span>
    <span class="n">band</span> <span class="o">=</span> <span class="n">raster</span><span class="o">.</span><span class="n">GetRasterBand</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">band</span><span class="o">.</span><span class="n">GetNoDataValue</span><span class="p">()</span>

<span class="k">def</span> <span class="nf">array2raster</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">,</span><span class="n">newRasterfn</span><span class="p">,</span><span class="n">array</span><span class="p">):</span>
    <span class="n">raster</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">Open</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">)</span>
    <span class="n">geotransform</span> <span class="o">=</span> <span class="n">raster</span><span class="o">.</span><span class="n">GetGeoTransform</span><span class="p">()</span>
    <span class="n">originX</span> <span class="o">=</span> <span class="n">geotransform</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
    <span class="n">originY</span> <span class="o">=</span> <span class="n">geotransform</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span>
    <span class="n">pixelWidth</span> <span class="o">=</span> <span class="n">geotransform</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
    <span class="n">pixelHeight</span> <span class="o">=</span> <span class="n">geotransform</span><span class="p">[</span><span class="mi">5</span><span class="p">]</span>
    <span class="n">cols</span> <span class="o">=</span> <span class="n">raster</span><span class="o">.</span><span class="n">RasterXSize</span>
    <span class="n">rows</span> <span class="o">=</span> <span class="n">raster</span><span class="o">.</span><span class="n">RasterYSize</span>

    <span class="n">driver</span> <span class="o">=</span> <span class="n">gdal</span><span class="o">.</span><span class="n">GetDriverByName</span><span class="p">(</span><span class="s">&#39;GTiff&#39;</span><span class="p">)</span>
    <span class="n">outRaster</span> <span class="o">=</span> <span class="n">driver</span><span class="o">.</span><span class="n">Create</span><span class="p">(</span><span class="n">newRasterfn</span><span class="p">,</span> <span class="n">cols</span><span class="p">,</span> <span class="n">rows</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="n">gdal</span><span class="o">.</span><span class="n">GDT_Float32</span><span class="p">)</span>
    <span class="n">outRaster</span><span class="o">.</span><span class="n">SetGeoTransform</span><span class="p">((</span><span class="n">originX</span><span class="p">,</span> <span class="n">pixelWidth</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="n">originY</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="n">pixelHeight</span><span class="p">))</span>
    <span class="n">outband</span> <span class="o">=</span> <span class="n">outRaster</span><span class="o">.</span><span class="n">GetRasterBand</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span>
    <span class="n">outband</span><span class="o">.</span><span class="n">WriteArray</span><span class="p">(</span><span class="n">array</span><span class="p">)</span>
    <span class="n">outRasterSRS</span> <span class="o">=</span> <span class="n">osr</span><span class="o">.</span><span class="n">SpatialReference</span><span class="p">()</span>
    <span class="n">outRasterSRS</span><span class="o">.</span><span class="n">ImportFromWkt</span><span class="p">(</span><span class="n">raster</span><span class="o">.</span><span class="n">GetProjectionRef</span><span class="p">())</span>
    <span class="n">outRaster</span><span class="o">.</span><span class="n">SetProjection</span><span class="p">(</span><span class="n">outRasterSRS</span><span class="o">.</span><span class="n">ExportToWkt</span><span class="p">())</span>
    <span class="n">outband</span><span class="o">.</span><span class="n">FlushCache</span><span class="p">()</span>


<span class="n">rasterfn</span> <span class="o">=</span> <span class="s">&#39;Slope.tif&#39;</span>
<span class="n">newValue</span> <span class="o">=</span> <span class="mi">0</span>
<span class="n">newRasterfn</span> <span class="o">=</span> <span class="s">&#39;SlopeNew.tif&#39;</span>

<span class="c"># Convert Raster to array</span>
<span class="n">rasterArray</span> <span class="o">=</span> <span class="n">raster2array</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">)</span>

<span class="c"># Get no data value of array</span>
<span class="n">noDataValue</span> <span class="o">=</span> <span class="n">getNoDataValue</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">)</span>

<span class="c"># Updata no data value in array with new value</span>
<span class="n">rasterArray</span><span class="p">[</span><span class="n">rasterArray</span> <span class="o">==</span> <span class="n">noDataValue</span><span class="p">]</span> <span class="o">=</span> <span class="n">newValue</span>

<span class="c"># Write updated array to new raster</span>
<span class="n">array2raster</span><span class="p">(</span><span class="n">rasterfn</span><span class="p">,</span><span class="n">newRasterfn</span><span class="p">,</span><span class="n">rasterArray</span><span class="p">)</span>
</pre></div>
</div>
</div>
</div>


          </div>
        </div>
      </div>
      <div class="sphinxsidebar" role="navigation" aria-label="main navigation">
        <div class="sphinxsidebarwrapper">
  <h3><a href="index.html">Table Of Contents</a></h3>
  <ul>
<li><a class="reference internal" href="#">Raster Layers</a><ul>
<li><a class="reference internal" href="#close-a-raster-dataset">Close a raster dataset</a></li>
<li><a class="reference internal" href="#get-raster-metadata">Get Raster Metadata</a></li>
<li><a class="reference internal" href="#get-raster-band">Get Raster Band</a></li>
<li><a class="reference internal" href="#loop-through-all-raster-bands">Loop Through All Raster Bands</a></li>
<li><a class="reference internal" href="#get-raster-band-information">Get Raster Band Information</a></li>
<li><a class="reference internal" href="#polygonize-a-raster-band">Polygonize a Raster Band</a></li>
<li><a class="reference internal" href="#convert-an-ogr-file-to-a-raster">Convert an OGR File to a Raster</a></li>
<li><a class="reference internal" href="#clip-a-geotiff-with-shapefile">Clip a GeoTiff with Shapefile</a></li>
<li><a class="reference internal" href="#calculate-zonal-statistics">Calculate zonal statistics</a></li>
<li><a class="reference internal" href="#raster-to-vector-line">Raster to vector line</a></li>
<li><a class="reference internal" href="#create-raster-from-array">Create raster from array</a></li>
<li><a class="reference internal" href="#create-least-cost-path">Create least cost path</a></li>
<li><a class="reference internal" href="#replace-no-data-value-of-raster-with-new-value">Replace No Data Value of Raster with new value</a></li>
</ul>
</li>
</ul>

  <h4>Previous topic</h4>
  <p class="topless"><a href="vector_layers.html"
                        title="previous chapter">Vector Layers</a></p>
  <h4>Next topic</h4>
  <p class="topless"><a href="projection.html"
                        title="next chapter">Projection</a></p>
  <div role="note" aria-label="source link">
    <h3>This Page</h3>
    <ul class="this-page-menu">
      <li><a href="_sources/raster_layers.txt"
            rel="nofollow">Show Source</a></li>
    </ul>
   </div>
<div id="searchbox" style="display: none" role="search">
  <h3>Quick search</h3>
    <form class="search" action="search.html" method="get">
      <input type="text" name="q" />
      <input type="submit" value="Go" />
      <input type="hidden" name="check_keywords" value="yes" />
      <input type="hidden" name="area" value="default" />
    </form>
    <p class="searchtip" style="font-size: 90%">
    Enter search terms or a module, class or function name.
    </p>
</div>
<script type="text/javascript">$('#searchbox').show(0);</script>
        </div>
      </div>
      <div class="clearer"></div>
    </div>
    <div class="related" role="navigation" aria-label="related navigation">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="genindex.html" title="General Index"
             >index</a></li>
        <li class="right" >
          <a href="projection.html" title="Projection"
             >next</a> |</li>
        <li class="right" >
          <a href="vector_layers.html" title="Vector Layers"
             >previous</a> |</li>
        <li class="nav-item nav-item-0"><a href="index.html">Python GDAL/OGR Cookbook 1.0 documentation</a> &raquo;</li> 
      </ul>
    </div>
    <div class="footer" role="contentinfo">
        &copy; Copyright 2013, Jared Erickson, Cort Daniel, Michael Payne.
      Created using <a href="http://sphinx-doc.org/">Sphinx</a> 1.3.1.
    </div>
  </body>
</html>